Loading Data

cod_data = read_csv("./data/NCHS_-_Potentially_Excess_Deaths_from_the_Five_Leading_Causes_of_Death.csv") %>%
  clean_names() %>%
  na.omit() %>%
  filter(!(state == "United States")) %>%
  separate(., percent_potentially_excess_deaths, into = c("percent_excess_death"), sep = "%") %>% 
  mutate(percent_excess_death = as.numeric(percent_excess_death), mortality = observed_deaths/population * 10000, mortality = as.numeric(mortality)) %>% 
  select(year, age_range, cause_of_death, state, locality, observed_deaths, population, expected_deaths, potentially_excess_deaths, percent_excess_death, mortality, hhs_region)
## Warning: Too many values at 191748 locations: 1, 2, 3, 4, 5, 6, 7, 8, 9,
## 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...

##columns removed

 #"state_fips_code"                      "benchmark"   "potentially_excess_deaths" "percent_excess_death"      "mortality"   

• the year variable contains data collected from 2005-2015.
• filtered out United States in the state variable.
• Added a variable mortality which is calculated by observed_deaths/population * 10000. This variable indicates the number of deathes observed in every 10000 people in the three geographic regions: Metropolitan, Nonmetropolitan and All.

region_cod_data = cod_data %>%
  select(state, locality, hhs_region, percent_excess_death) %>% 
  group_by(state,locality, hhs_region) %>% 
  summarise(mean_ped = mean(percent_excess_death)) %>% 
  dplyr::filter(!(state == "District of\nColumbia")) %>% 
  mutate(hhs_region = as.character(hhs_region))


region_lm = lm(region_cod_data$mean_ped~region_cod_data$hhs_region)
rl1 = broom::tidy(region_lm)
kable(rl1)

region_locality_lm =lm(region_cod_data$mean_ped ~ region_cod_data$hhs_region+region_cod_data$locality)
rl2 = broom::tidy(region_locality_lm)
kable(rl2)

Interpretation: U.S. Department of Health and Human Services public health regions (1 through 10) are used as a categorical variable in the above regressions. Specific region classification is shown in the figure below.

[attach image here].

As shown in the above summary tables, with respect to region 1, only region 2 have negative estimated coefficient, indicating less mean percentage of excess death . Specifically, comparing with region 1, region 2 have 1.67% less mean percentage of excess death on average. From region 3 to region 10, the mean percentage of excess death is higher comparing with region 1. Similiar results yield from regression adjusted for Locality. Adjusting locality, region 4 and 6 have top two highest increase in mean percetage of excess death with respect to region 1. Adjusting for different regions, similar result yield that, on average, metropolitan have 2% less than overall mean percetage of excess deaths and nonmetropolitan have 7% more than overall mean percetage of excess death.

As shown in the above summary tables, with respect to region 1, only region 2 have negative estimated coefficient, indicating less mean percentage of excess death . Specifically, comparing with region 1, region 2 have 1.67% less mean percentage of excess death on average. From region 3 to region 10, the mean percentage of excess death is higher comparing with region 1. Similiar results yield from regression adjusted for Locality. Adjusting locality, region 4 and 6 have top two highest increase in mean percetage of excess death with respect to region 1. Adjusting for different regions, similar result yield that, on average, metropolitan have 2% less than overall mean percetage of excess deaths and nonmetropolitan have 7% more than overall mean percetage of excess death.

Locality vs. Mortality

cod_data %>%
  mutate(year = as.factor(year))
## # A tibble: 191,748 x 12
##      year age_range cause_of_death   state        locality observed_deaths
##    <fctr>     <chr>          <chr>   <chr>           <chr>           <int>
##  1   2005      0-49         Cancer Alabama             All             756
##  2   2005      0-49         Cancer Alabama    Metropolitan             556
##  3   2005      0-49         Cancer Alabama Nonmetropolitan             200
##  4   2005      0-49         Cancer Alabama             All             756
##  5   2005      0-49         Cancer Alabama    Metropolitan             556
##  6   2005      0-49         Cancer Alabama Nonmetropolitan             200
##  7   2005      0-49         Cancer Alabama             All             756
##  8   2005      0-49         Cancer Alabama    Metropolitan             556
##  9   2005      0-49         Cancer Alabama Nonmetropolitan             200
## 10   2005      0-54         Cancer Alabama             All            1346
## # ... with 191,738 more rows, and 6 more variables: population <int>,
## #   expected_deaths <int>, potentially_excess_deaths <int>,
## #   percent_excess_death <dbl>, mortality <dbl>, hhs_region <int>
p <- cod_data %>%
  plot_ly(
    x = ~expected_deaths, 
    y = ~observed_deaths, 
    size = ~population, 
    color = ~cause_of_death, 
    frame = ~hhs_region, 
    text = ~state, 
    hoverinfo = "text",
    type = 'scatter',
    mode = 'markers'
  ) %>% 
  layout(title = "Change of Standardized Mortality Ratio in National Public Health Regions")

p

cod_data %>%
  group_by(cause_of_death) %>% 
  ggplot(aes(x = locality, y = mortality, fill = cause_of_death)) +
  geom_col() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title="Locality vs. Mortality") +
  labs(x="locality", y="mortality") 

cod_data %>%
  group_by(cause_of_death) %>% 
  ggplot(aes(x = locality, y = mortality)) + geom_boxplot(aes(color = cause_of_death), na.rm = T) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title="Locality vs. Mortality") +
  labs(x="locality", y="mortality") 

mortality_lm = lm(mortality ~locality, data = cod_data)
summary(mortality_lm)
## 
## Call:
## lm(formula = mortality ~ locality, data = cod_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1454 -2.9547 -1.2983  0.9988 19.8959 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              4.10331    0.01672 245.474   <2e-16 ***
## localityMetropolitan    -0.22576    0.02372  -9.519   <2e-16 ***
## localityNonmetropolitan  1.14451    0.02438  46.942   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.317 on 191745 degrees of freedom
## Multiple R-squared:  0.01821,    Adjusted R-squared:  0.0182 
<<<<<<< HEAD
## F-statistic:  1778 on 2 and 191745 DF,  p-value: < 2.2e-16
rl3 = broom::tidy(mortality_lm)
kable(rl3)

Mean Percent Excess Death vs. Locality

gp_cod_data = cod_data %>%
  group_by(year, locality) %>% 
  summarise(mean_ped = mean(percent_excess_death)) 

renderPlotly({
gp_cod_data_Metro = gp_cod_data %>%
  group_by(year) %>%
  filter(locality == "Metropolitan")
gp_cod_data_nonMetro = gp_cod_data %>%
  group_by(year) %>%
  filter(locality == "Nonmetropolitan")
gp_cod = cbind(gp_cod_data_Metro, gp_cod_data_nonMetro)

p = gp_cod %>%
plot_ly( x = ~year, y = ~mean_ped, name = 'Metropolitan', type = 'scatter', mode = 'markers') %>%
  add_trace( y = ~mean_ped1, name = 'Nonmetropolitan', mode = 'markers') %>% 
  layout(yaxis = list(title = 'mean percent excess death (Metropolitan/Nonmetropolitan)'))
p
})
## function (...) 
## {
##     if (length(outputArgs) != 0 && !hasExecuted$get()) {
##         warning("Unused argument: outputArgs. The argument outputArgs is only ", 
##             "meant to be used when embedding snippets of Shiny code in an ", 
##             "R Markdown code chunk (using runtime: shiny). When running a ", 
##             "full Shiny app, please set the output arguments directly in ", 
##             "the corresponding output function of your UI code.")
##         hasExecuted$set(TRUE)
##     }
##     if (is.null(formals(origRenderFunc))) 
##         origRenderFunc()
##     else origRenderFunc(...)
## }
<<<<<<< HEAD
## <environment: 0x7faa0ae5aae0>
=======
## <environment: 0x7ffdcf4876d0>
>>>>>>> 667868adab929cae16d0a55c63183dce21b1a2c5
## attr(,"class")
## [1] "shiny.render.function" "function"             
## attr(,"outputFunc")
## function (outputId, width = "100%", height = "400px", inline = FALSE) 
## {
##     htmlwidgets::shinyWidgetOutput(outputId = outputId, name = "plotly", 
##         width = width, height = height, inline = inline, package = "plotly")
## }
## <environment: namespace:plotly>
## attr(,"outputArgs")
## list()
## attr(,"hasExecuted")
## <Mutable>
##   Public:
##     clone: function (deep = FALSE) 
##     get: function () 
##     set: function (value) 
##   Private:
##     value: FALSE

gp_locality_lm = lm(mean_ped ~locality , data = gp_cod_data)
summary(gp_locality_lm)
## 
## Call:
## lm(formula = mean_ped ~ locality, data = gp_cod_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5667 -1.0325 -0.3470  0.7922  3.2937 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              33.8866     0.4071  83.238  < 2e-16 ***
## localityMetropolitan     -2.0599     0.5757  -3.578   0.0012 ** 
## localityNonmetropolitan   8.0460     0.5757  13.975 1.13e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.35 on 30 degrees of freedom
## Multiple R-squared:  0.9198, Adjusted R-squared:  0.9145 
## F-statistic: 172.1 on 2 and 30 DF,  p-value: < 2.2e-16
rl4 = broom::tidy(gp_locality_lm)
kable(rl4)

Mean Percent Excess Death vs. Causes of Death

How should we help nonmetropolitan areas improve public health?

cod_datayr2 = read_csv("./data/NCHS_-_Potentially_Excess_Deaths_from_the_Five_Leading_Causes_of_Death.csv") %>%
  clean_names() %>%
  na.omit() %>%
  filter(!(state == "United States")) %>%
  separate(., percent_potentially_excess_deaths, into = c("percent_excess_death"), sep = "%") %>%
  mutate(percent_excess_death = as.numeric(percent_excess_death), mortality = observed_deaths/population * 10000, mortality = as.numeric(mortality)) %>%
  select(year, cause_of_death, state, benchmark, locality, observed_deaths, population, expected_deaths, potentially_excess_deaths, percent_excess_death, mortality)
## Warning: Too many values at 191748 locations: 1, 2, 3, 4, 5, 6, 7, 8, 9,
## 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...

cod_datayr2 %>%
  filter(benchmark == "2010 Fixed")%>%
  group_by(cause_of_death,year, locality)%>%
  mutate(year_percent_death=mean(percent_excess_death))%>%
  distinct(year, cause_of_death, locality, .keep_all = TRUE)%>%
  ungroup(year, locality)%>%
ggplot(aes(x = year, y = year_percent_death, color = cause_of_death))+
  geom_point(size = 1)+
  geom_line(size = 1) +
  facet_grid(.~locality)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  scale_x_continuous( breaks=seq(2005,2015,1))+
  labs(title = "Average yearly excess death percentage for different locality")

interpretation on linear regression: As can be seen in the plot, if we stratify by cause, we will find some increasing pattern in subgroups(unintentional injury) in certain period of time(2010-2015). In nonmetropolitan regions, unintentional injury and chronic lower respiratory disease have on average the highest excess death percentage in nonmetropolitan areas and needs more attention from public health experts.

cod_datayrlm = cod_datayr2 %>%
  filter(locality == "Nonmetropolitan",
         benchmark == "2010 Fixed")
lmpop = lm(data=cod_datayrlm, percent_excess_death~year+cause_of_death)
summary(lmpop)
## 
## Call:
## lm(formula = percent_excess_death ~ year + cause_of_death, data = cod_datayrlm)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -55.456  -9.828   0.576  10.404  39.489 
## 
## Coefficients:
##                                                  Estimate Std. Error
## (Intercept)                                     769.06335   67.48938
## year                                             -0.36941    0.03358
## cause_of_deathChronic Lower Respiratory Disease  27.06490    0.33901
## cause_of_deathHeart Disease                      13.90257    0.32844
## cause_of_deathStroke                             12.91441    0.33784
## cause_of_deathUnintentional Injury               31.18470    0.32812
##                                                 t value Pr(>|t|)    
## (Intercept)                                       11.39   <2e-16 ***
## year                                             -11.00   <2e-16 ***
## cause_of_deathChronic Lower Respiratory Disease   79.83   <2e-16 ***
## cause_of_deathHeart Disease                       42.33   <2e-16 ***
## cause_of_deathStroke                              38.23   <2e-16 ***
## cause_of_deathUnintentional Injury                95.04   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.92 on 19718 degrees of freedom
## Multiple R-squared:  0.365,  Adjusted R-squared:  0.3649 
## F-statistic:  2267 on 5 and 19718 DF,  p-value: < 2.2e-16
rl5 = broom::tidy(lmpop)
kable(rl5)

From the linear model based on 2010 fixed benchmark, we can tell that controlling for cause of death, the excess death percentage is expected to decrease by 0.37 each year. Overall, the excess death percentage shows a decreasing pattern with time. Besides, in the same year, compared with cancer, other 4 causes have significantly higher excess death percentage. Herein unintentional injury and chronic lower respiratory disease have on average the highest excess death percentage.

Mean Percentage of Excess Death vs. Five Leading Cause of Death

For each of the five cause of death, what is the mean percent of excess death? It is different across locality. As shown in the boxplots for all three regions, the rank of the five cause of death remain the same for different localities. For all five causes, the mean percent of death is lower for Metropolitan and higher for Nonmetropolotan.

For each of the five cause of death, what is the mean percent of excess death? It is different across locality.

<<<<<<< HEAD As shown in the boxplots for all three regions, the rank of the five cause of death remain the same for different localities. For all five causes, the mean percent of death is lower for Metropolitan and higher for Nonmetropolotan.

cod_data %>%
  select(state, locality, percent_excess_death, hhs_region, cause_of_death) %>%
  filter(locality != "All") %>% 
  mutate(hhs_region = as.factor(hhs_region)) %>% 
  group_by(cause_of_death, locality) %>% 
  summarise(mean_ped = mean(percent_excess_death)) %>% 
  ungroup() %>% 
  ggplot(aes(x = cause_of_death,  y = mean_ped, fill = cause_of_death)) +
  geom_col() +
  facet_grid(~locality) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 

<<<<<<< HEAD

For the above graph, if we seperate the effect of the gap between urban and rural area into different HHS region(Figure “above”) as mentioned, the following chart gives detailed information.

cod_data %>%
  select(state, locality, percent_excess_death, hhs_region, cause_of_death) %>%
  filter(locality != "All") %>% 
  mutate(hhs_region = as.factor(hhs_region)) %>% 
  group_by(cause_of_death, locality,hhs_region) %>% 
  summarise(mean_ped = mean(percent_excess_death)) %>% 
  ggplot(aes(x = hhs_region,  y = mean_ped, fill = hhs_region)) +
  geom_col() +
  facet_grid(~cause_of_death) +
  theme(axis.text.x = element_text(angle = 15)) +
    labs(
      y = "Mean %Excess Death",
      x = "HHS Regions",
      fill = "HHS Region"
    ) 

In order to take the effect of Localtiy into account, Flexdashboard with Shiny will be helpful. the link to the Shiny page is Here - “Comparison of Difference in Mean %Excess Death of Five Leading Causes of Death in each HHS regions between Metropolitan and Nonmetropolitan Area”.

Regress percent of excess death on cause of death, which is a categorical variable. with respect to cancer, the least of the five causes, all 4 have additional percent of death as evaluated by the estimates. The results give statisically significant estimates .

interpretation on linear regression: yz3306:Regress percent of excess death on cause of death, which is a categorical variable. with respect to cancer, the least of the five causes, all 4 have additional percent of death as evaluated by the estimates. The results give statisically significant estimates.

cod_lm_yz1 = lm(percent_excess_death ~cause_of_death, data = cod_data)
rl6 = broom::tidy(cod_lm_yz1)
kable(rl6)